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Abstract 

In many applications, such as those arising from the field of cellular networks, it is often desired to 
determine the interaction (graph) structure of a set of differential equations, using as data measured 
sensitivities. This note proposes an approach to this problem. 

1 Introduction 

Suppose given a system of differential equations 

x\ = fi(xi,...,x n ,pi,...,Pm) 

X2 = h{x\, ■ ■ ■ ,X n ,Pl, ■ ■ ■ ,Pm) 

■. w 

X71 — fn \X\ 1 • • • 5 X n , p\ , . . . , Pro) , 

where the vector of state variables x(t) = (xi(t), . . . ,x n (t)) evolves in some open subset X C M™ and 
the vector of parameters p = (pi, . . . ,p m ) can be chosen from an open subset V C R m , for some n and 
7ti. In cellular networks, for instance, the state variables Xi might represent the concentrations of certain 
proteins, mRNA, etc., at different time instants, and the parameters^ might represent the concentration 
levels of certain enzymes which are maintained at a constant value during a particular experiment (see, 
e.g., 0,1). 

In many fields of application, it is often the case that the equations defining the system (that is, the 
form of the functions fi describing the vector field) are unknown, even in general form, but one wishes 
nonetheless to determine the interaction graph of a system (TJ) , that is to say, to know which variables 
directly influence which variables, as well as the relative strengths of these interactions (|3j). A more 
limited goal might be to determine if a certain feedback loop is present in the system (cf. [4 ). 

To help in this task, experimental data is usually available, measuring solutions of system (|l|) for 
various initial states and parameter settings. A special case, the one treated in this note, is that in which 
experimental data provides us with the location of steady states associated to parameter values near a 
given set of parameters p. We show in this note how, starting from such data, and under an assumption 
that seems natural, it is indeed possible to determine this interaction graph and relative strengths of 
interactions. We then extend the method to non-steady state measurements. A more global analysis of 
the problem is also possible, and will be presented in a follow-up report. 
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2 Problem Formulation 



The steady states of system (0) associated to a given parameter vector p = (pi, . . . ,p m ) are the solutions 
x = (x\, . . . , x n ) of the set of n simultaneous algebraic equations 

fl(xi,...,X n ,pi,...,p m ) = 

f 2 {x 1 ,...,x n ,p 1 ,...,p m ) = 

: ( 2 ) 

f n (xi,...,X n ,p 1 ,...,p m ) = 0. 

We will assume that there is a function £ : V — » X which assigns, to each parameter vector p G V, a 
steady state of (jl]), that is to say, one has that fi(£(p),p) — for all i = 1, . . . ,n and all p G P. 
We suppose that the functions fi as well as £ are continuously differentiable. A particular parameter 
vector p is also given, and the problem will be formulated in terms of the behavior of steady states near 
x = f(p). 

Experimental Data 

It will be assumed that an n x m matrix E = (akj) is given, representing the "sensitivities" 

^ = ||(P) (3) 

for each k — 1, . . . ,n and each j — 1, . . . , m. This matrix of partial derivatives may be estimated 
numerically from the values of £ (p) on a neighborhood of the chosen parameter p. 

Desired Information 

Consider the n x n matrix A = (a^) defined by: 
for every i, j = 1, . . . , , n. 

Ideally, one would want to find the matrix A, since this matrix completely describes the influence of 
each variable Xj upon the rate of change of each other variable Xi. Unfortunately, such an objective is 
impossible to achieve from the local steady-state data E, or even from the knowledge of the complete 
(global) mapping £. This is because the same mapping £ also solves the set of equations 



X 1 f 1 (x i ,...,x n ,p 1 ,...,p m ) = 

^2f2(xi,...,X n ,pi,...,p m ) = 
fn j • • • ) %n i Pi •>••••> Pm ) — 5 



(4) 



for any constants A^, but, on the other hand, multiplication of fi by Xi results in a^- = XiCiij. In other 
words, the best that one could hope is for the data E to determine the rows 

A t = (an, . . . , din) , i = 1, ... ,n 

of A only up to scalar multiples. 

Thus, a more realistic objective is to attempt to identify the rows Ai up to a scalar multiple only. 
For example, if we assume that an ^ for each i (a realistic assumption when stable systems are being 
interconnected), this amounts to finding the ratios aij/au for each i =/= j. 
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Assumptions 



We will make two assumptions which will suffice for us to solve the problem of determining the rows Aj 
of A up to scalar multiples. The first assumption is a strong but reasonable structural one, while the 
second represents a weak algebraic nondegeneracy condition. 

We will suppose known, for each i S {1, . . . , n}, a subset Si of the index set {1, . . . , to} so that the 
following property holds: 

(VjeSi) M(s,p) = Q (5) 
d Pj 

which is in turn implied by the structural condition: 

(y j G Si) fi does not depend upon pj . 

This prior information about the system structure is far less restrictive than it might appear at first 
sight. Indeed, it is usually the case that "compartmcntal" information is available, for instance telling us 
that the concentration of a certain enzyme has no direct influence on an unrelated biochemical reaction, 
that an extracellular signaling molecule does not affect directly a cytoplasmatic reaction, and so forth. 

The second assumption is as follows. For each j G {!,..., to}, we introduce the vector 



Sj = col (ctij, . . . , a n j) 



HP) 



(6) 



representing the jth column of the matrix E, and we consider, for each i G {1, . . . , n} the linear subspace 
Hi of K n spanned by the vectors 

The assumption is: 

dim Hi > n — 1 V i = 1, . . . , n . 

Note that this amounts to saying that the dimension of Hi is either n or n — 1. (Generically, we may 
expect this dimension to be n — 1, because the orthogonality relation to be shown below is the only 
algebraic constraint.) 



3 Solution 

With the above assumptions, the problem that we posed can be solved as follows. We fix any index 
i G {1, . . . , n}, and take partial derivatives in the equation fi(x(p),p) = with respect to the variable 
Pj, for each index j in the set Si, and evaluate at x = x and p — p: 




where the second term vanishes by assumption (J^). 

Since this happens for every j G Si, we conclude that the vector Ai is orthogonal to Hi. As Hi has 
dimension n or n — 1, this determines Aj up to a scalar multiple, which is what we wanted to prove. 

Of course, it is trivial to actually compute Ai from the data. If Hi has dimension n, then Ai = 0; 
this is a degenerate case. When the dimension is n — 1, one simply picks a basis {E^ , . . . , Sj n _ 1 } of Hi, 
and any vector E linearly independent from the elements of this basis (a randomly chosen vector has 
this property), and then solves (for example) the nonsingular set of equations Aj • Eg = 1, Aj ■ E Jf = 0, 
£ = I, . . . ,n — 1, to find a nonzero Aj (all possible Aj are scalar multiples of this one). Alternatively, 
provided that one knows that an 0, one may simply normalize to an = —1 and then determine the 



remaining entries of Ai by solving a linear set of n — 1 equations. Observe also, that if it is known a 
priori that certain entries vanish, then one may redefine the space Hi to be spanned only by the 
vectors listing the appropriate components of the sensitivities d^/dp^s, and a potentially much smaller 
number of parameter perturbations may be required. 



4 Modular Approach 

It is also possible to apply our techniques in a "modular" context, in which only the derivatives dfi/dxj 
with respect to communicating intermediaries are calculated (Q). Let us briefly explain this. 

We assume that the entire network consists of an interconnection of n subsystems or "modules" , each 
of which is described by a set of differential equations such as: 



Xj 


= 9oj(yi,j,- 


-,ye,j,xi, ■ ■ 


■ 7 X n ,Pl,.. 


■ ■ ,Pm) 


Vl,3 


= gx,j(yx,j,.. 


■ ,yt,j,xi, . . 


■ j x n , pi , . . 


■ ■ ,Pm) 


V2,j 


= 92,j{yi,j,-- 


■ ,yt,j,xi, . . 


■ j x n ,pi,.. 


■ ■ ,Pni) 



ytjd = 9t j ,j(yi,j>---,y£,j> x i>---> x n,pi,---,Pm), 

where the variables Xj represent "communicating" or "connecting" intermediaries of module j that 
transmit information to other modules, whereas the variables yi.j, ■ ■ ■ ,Vtj represent chemical species 
that interact within module j. The integer £j, j = 1, . . . , n is in general different for each of the n 
modules and represents the number of chemical species in the jth module. 

We will assume that, for each fixed module, the Jacobian of (<7i, . . . , gi, ( j ) with respect to y\, . . . , yi . , 
evaluated at the steady state corresponding to p (assumed to exist, as before) is nonsingular. The 
Implicit Mapping Theorem then implies that one may, in a neighborhood of this steady state, solve 

90, j(yi,j,---,yi,j,xi,...,x n ,Pu--.,Pm) = o 

91, j(yi,j,---,ye,j,xx,...,x n ,pi,...,p m ) = o 

92, j(yi,j,---,yi,j,xx,...,x n ,pi,...,p m ) = o 

9e j ,j(yi,j,---,ye,j,xu---,Xn,Pi,---,Pm) = 

for the variables Xj, y%, . . . , yi j as a function of x\, . . . , x ni p\ 1 . . . ,p m . One concludes that, around this 
steady state corresponding to p, the functions Xj satisfy implicit equations of the form 

Xj hj (xi , . . . , x n , pi , . . . , p ra ) 

which we can rewrite in the form (Q), using fj(x,p) = Xj — hj(x,p). The analysis then proceeds as 
before. The generalization to the case of more than one communicating intermediate in a module, 
namely a vector (xj t i, . . . , is obvious. 



5 Avoiding Derivatives 

The technique that was described assumes that we know the sensitivity matrix S, which is obtained 
by evaluating the partial derivatives d^i/dpj at the particular parameter value p. Ordinarily, these 
derivatives would be estimated by finite differences. For instance, suppose that one measures x, = 
as well as £(p + djp), where 

djp = col (0, . . . , 0, dpj , 0, . . . , 0) 
(entry in jth position) and we view dpj as a "small" perturbation of the jth parameter. Denoting 

d 3 Xi := £i(p+ djp) - Xi (9) 
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obviously one may estimate S using the following approximation: 



-(P) 



dpj dpj 

In order to calculate this ratio, both djXi and dpj must be known. 

However, in certain experimental situations it may well be impossible to estimate the values of dpj. 
This might appear to be contradictory, since we are assuming that we perform experiments which change 
the values of p. But one can easily envision an experimental setup in which a certain external variable 
(in a cell biology situation, for instance, a growth factor) is known to influence a certain parameter pj. 
Varying this external variable therefore produces a perturbation in pj, and hence an appropriate djX 
which is measured, but dpj itself may be hard to measure. 

It is rather surprising that we can still achieve our goal of estimating the rows of A (up to scalar 
multiples) even in the absence of information about the dp^sl To see intuitively why this is plausible, 
consider the following argument. Let us say that we have just a scalar parameter p and a scalar function 
f(xi,X2) so that f(£,i(p),£,2{p)) = 0, and that x — £(p) — (0,0). In a neighborhood of p = p, we may 
assume that / is linear, so we have a linear relation (with unknown coefficients) 

Ol6(p) +02&(P) = 0. 

The method discussed so far would take derivatives at p = p: 

oi-t-(p) + «2-j-(p) = 
dp dp 

and thus (assuming that the derivative is not zero), we know that the row (01,02) must be a multiple 
of {—{d£ t 2/dp)(p),{d£,i/dp)(p)). A completely different argument (analogous to using a two-point as 
opposed to a slope-point formula in order to find the equation of a line) would simply take the original 
equation <Zi£i(f>) + 02^2 (p) = (valid only p w p, since this was an approximation of /) and say that 
the row (01,02) must be a multiple of {—^{p), £i(p))> for any fixed p « p. There is no inconsistency 
between the two estimates, since they only differ (approximately) by multiplication by the scalar dp: 

(-€a(p),6(p)) « (-(^/dp)(p),(d^/dp)(p))dp 
and we only care about scalar multiples. Let us now say this in general. 

Since /i(£(p + djp),p + djp) — f(x,p) = — = 0, we have, taking a Taylor expansion, that 

d " ^ Q j. 

-^fi(€(p),P) _JPj + o(dpj) = — (p)dj?j + o{d Pj ) = 0. (10) 



— iln.. ' dpj 



whenever j £ Si. Substituting 



djx k = £ k (p + djp) - x k = -^—{p)dp.j + o(dpj) (II) 

Opj 



into (JToj) , we conclude that 



" d fi 

} -—^-(x,p)djXk — o(dpj) provided that j G Si . 

, oxk 

Since dpj « and djXk = 0(dpj), this is an approximate orthogonality relation, and we now make the 
approximation: 

} -^—?-{x,p)djXk — provided that j £ Si . (12) 

In conclusion, and introducing the matrix Y = (jkj) = {djXk) instead of S, we have that A4 ■ Tj = for 
all j G Si, where Yj = col(7i,-, . . . ,7 n j)- Now we consider, for each i G {1, . . . the linear subspace 
Ki of W 1 spanned by the vectors {r^ \j G Si} and assume that dim Ki > n — 1 for all i. We conclude 
that the vector A4 is orthogonal to Ki, and this once again determines Ai up to a scalar multiple. 
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6 Non- Steady State Analysis 



Let us sketch here how one might extend our methodology to use non-steady state data. In general, we 
denote by £(t, x°,p) the solution of (Q) with initial condition a; , at time t and using parameters p. Let 
us suppose that we can measure the sensitivities (|^) at some specific point in time, and for some specific 
solution £(t, x°,p): 

a kj = ^(i,x°,p) (13) 

for each k = 1, . . . ,n and each j = 1, . . . , m, and we let £ = (c/cj)- We also need now the mixed second 
derivatives: 

^^dp-k^ ^ 

and instead of Ej = col (cry, . . . , cr n j) as m @j we consider for each i G {1, . . . , n} and j S {1, . . . , m} 
the vector 

Eij = co^ryy,^, . . .,<r„j). (14) 

We define, for each i 6 {1, . . . , n}, Hi as the linear subspace of spanned by the vectors {S^- | j £ Si}. 
We let now = x°,p). Fixing any index i £ {1, . . . ,n}, we take partial derivatives on both sides 

of the differential equation 

^i(t,x°,p) = fMt,x a ,p),p) 

with respect to the variable pj, for each index j in the set Si, and evaluate at x — x, t — i, and p = p, 
to obtain: 



Vij = -S—h{x,p) = y^a ik (7kj 



from which we conclude that [—1, Ai] ■ Yj } ■ — whenever j 6 Si, and hence that [—1, Ai] is orthogonal to 
Hi. With appropriate genericity conditions, this orthogonality, perhaps in conjunction with conditions 
at other times t or points p, will restrict the possible vectors Ai and more generally the interaction graph. 
(For example, if dim Hi — n, then we have a unique solution.) Derivatives with respect to parameter 
values can be replaces by differences, just as in the steady state case. We will discuss this further in a 
future contribution. 
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